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Abstract: Propagation of electromagnetic waves over irregular, inhomogeneous terrain is 
solved by a finite difference scheme. The method is fast and requires considerably less 
memory than the integral equation methods. The method requires a storage space of 
order O(N) and an execution time of order 0(N 2 ). Fields generated by a TE Z line source 
are represented in an integral form in terms of the field over a flat, constant impedance 
plane, and the field scattered by the terrain irregularities and inhmogeneities. Accurate 
expressions are provided for the incident field and the Green’s function, whose evaluation 
is otherwise accomplished by the rather time-consuming Sommerfeld’s integrals. Measured 
equation of invariance is used to terminate the mesh. The sparse matrix generated by the 
method is inverted by the Ricatti transform. Numerical results are presented for the ground 
wave as well as the sky wave. Comparison is made for known geometries to establish the 
validity and limitations of the method. 



I. INTRODUCTION 



Electromagnetic wave propagation over hilly terrain is important not only in point-to- 
point communication over land, but also in ground-to-air communication. Of late, it has 
assumed importance in outdoor propagation in the context of personal communications 
network design. An exact analytical solution of the problem for general terrain is not pos- 
sible and one often resorts to an approximate or a numerical approach. In a previous paper 
[1], we developed a numerical model for propagation predictions over inhomogeneous, ir- 
regular terrain using the magnetic field integral equation. Although the method could 
include all effects of wave propagation such as reflection, diffraction, surface wave excita- 
tion, and backscattering, a principal limitation of the method was the requirement of large 
computer resources (CPU time and memory), particularly for electrically large terrain ir- 
regularities. For instance, if the integral equation is solved numerically by the method of 
moments [2], the matrix fill time would be of order 0(A r2 ) and the inversion time of order 
0(jV 3 ), where N is the total number of unknowns. As the matrix generated is dense, the 
memory requirements would be of order 0(N 2 ). The method is attractive for small terrain 
irregularities, but is computationally prohibitive for large terrain irregularities [1]. 

In this paper we present a computationally efficient model of wave propagation based 
on finite differences. It may be noted that this method has no semblance to the one 
proposed in [3], where one proceeds with the parabolic equation approximation of the 
Helmholtz equation. Use of parabolic equation approximation precludes backscattering, 
which is sometimes important. In contrast, we apply finite differences directly to the 
Helmholtz equation without introducing any dubious approximations. Even though one 
has to handle a larger matrix when dealing with finite differences in contrast to boundary 
methods such as integral equation methods, the resulting matrix is sparse, and often faster 
to invert than the dense matrix generated via the latter. The matrix fill time is still of 
order 0(iV 2 ), but the inversion time is reduced to a lower order of 0(N). Thus, the real 
advantage of a finite difference scheme is felt for large problems where the run time is 
dominated by the inversion time. Substantial savings in memory are also be accomplished 
as the matrix is sparse. The memory required is only of order 0(N). The method can 
potentially solve larger problems than possible with the method of moments. The principal 
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difficulty associated with the use of finite differences is the treatment of mesh truncation 
when applied to open type of problems such as encountered in propagation, antennas, and 
scattering. We will use a method known as the Measured Equation of Invariance (MEI) 
originally proposed by Mei, et al, [4] to treat the mesh truncation. 

The term scatterer will be used generically to denote an obstacle in an open envi- 
ronment. Local radiation conditions such as those proposed in [5] and [6] are good when 
the truncating boundary is far removed so that it is in the far-zone field. However, this 
results in a larger matrix size with obvious implications to computational time and mem- 
ory. Global type radiation conditions [7] can be used on a tighter terminating boundary 
to result in a smaller computational domain. However, this will destroy the sparsity of the 
matrix generated by finite differences and defeats the whole purpose of employing it in the 
first place. What is needed is a local boundary condition of the type in [5], but applicable 
very close to the scatterer. Although it is far more complex to find near field radiation 
conditions than it is to find far field ones, it is partially offset by the fact that boundary 
conditions on a continuous spatial domain are not needed when finite methods are used. 
The MEI method enables one to generate the near-field conditions over a discrete domain. 

In this paper we deal with only two-dimensional sources and fields, and one dimen- 
sional terrain characteristics. Accordingly, the terrain properties, the sources, and the 
corresponding fields are all invariant with respect to the longitudinal variable 2 . It is as- 
sumed that the terrain is characterized by its local impedance and height over a reference 
plane, both of which may vary from point to point. In section Ha, we present a finite dif- 
ference discretization of the two dimensional Helmholtz equation and present an overview 
of the present method. To realize the mesh termination conditions via the MEI method 
described in lid, an accurate representation of the near-zone scattered field is necessary. 
In section lib we give an integral representation of the scattered field of a TE Z line source 
over an irregular, impedance surface. The corresponding expressions for the incident field 
and the Green’s functions are presented in section lie. In section He, we describe the Ri- 
catti block- by-block elimination technique [8] of sparse matrix inversion. Finally in section 
III we present numerical results for both the sky wave and the ground wave and provide 
comparisons for test geometries. 
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II. FORMULATION 



In the present paper, we consider a two dimensional situation as shown in Fig. 1. The 
transmitting antenna is a transversely polarized electric line source located at ( x 0 ,y 0 ). 
Such a source will have its electric field confined to the transverse xy plane, and is the 
two-dimensional counterpart of the superposition of a vertical electric dipole (VED) and 
a horizontal electric dipole (HED) in three dimensions. The field due to the source can 
be classified as TE 2 , and all components could be expressed in terms of the z-component, 
H z , of the magnetic field. It is assumed that all distance variables are normalized with 
respect to the free-space wavenumber k 0 = 0Jy/li o e oi where u is the radian frequency of 
the wave, e Q is the permittivity, and // 0 the permeability of free-space. Accordingly, we set 
x := k Q x , y := k 0 y , etc, An e JLji time dependence is assumed and suppressed. For TE 2 
polarization, the impedance boundary condition [9] of the form h x E = rj 0 Ah x (n x H ), 
relating the electric field vector E to the magnetic field vector H leads to 



dH z 

dn 



= ]AH 2 , 



(i) 



where the unit normal h points out of the impedance surface, t] q = yj j.i 0 / e 0 is the intrinsic 
impedance of free space, and A is the normalized impedance of the surface. 



Ha . Overview of the Method: 

We let ip to denote the z-component of the scattered magnetic field due to a TE 2 line 
source over an inhomogeneous, irregular terrain. We assume that rp represents scattering 
only from the irregularites and inhomogeneities in a reference impedance plane. Thus, 
the scattered field is identically zero when the terrain is flat having an impedance equal 
to the reference impedance. The scatterer then consists of those portions of the terrain 
where (i) the impedance is different from the reference value, and (ii) the elevation is 
different from zero. The computational domain consists of a region in space bounded by 
the terrain at the bottom and a terminating (i.e., truncating) boundary at the top as 
shown in Fig. 2. It is assumed that the boundary of the scatterer is subdivided into A r — 1 
segments, thereby, generating N points on it. The terminating boundary is similarly 
partitioned into N — 1 segments. We generate a structured mesh in the computational 
domain by adding M interior layers between the object boundary and the terminating 



3 






boundary, each having N points. A total of M + 2 layers, each having N points are 
thus generated. Layer numbering is done in ascending order starting from the object 
boundary (layer number 0) and progressing towards the terminating boundary. Node 
numbering is done from left (number 1) to right (number N). We use the notation (x^, y 
m = 0, 1 . . . , M + 1, n = 1, . . . , N to denote the cartesian coordinates of the nth point on 
layer number m. Points on the layer immediately following the scatterer are assumed to 
lie on the local normal emanating from a point on the scatterer. This is to permit easy 
implementation of the impedance boundary condition of (1). Within the computational 
domain, the scattered field \j) satisfies the scalar Helmholtz equation 

( d 2 d 2 \ 

(a^ + ^r + * = 0 - (2) 

Since the terrain irregularities do not necessaily conform to any one standard coordinate 
system, the mesh is non-orthogonal, and traditional finite difference equations are not 
applicable. We will develop the required finite difference equations similar to the method 
outlined in [10]. Prompted by the presence of second order derivatives in the Helmholtz 
equation, we choose a five-point star mesh centered at 0(Xo,Vo) and sorrounded by four 
neighboring nodes with local coordinates (X*, Y^), k = 1 , ... ,4, as shown in Fig. 3b. The 
global indices of the nodes 0, 1, 2, 3, and 4 are (m, n), (m — 1, n), (m + 1, n), (m, n — 1), and 
(m,n + 1) respectively. Using the notation in the local coordinates that V’fc = 
we assume a finite difference equation over the mesh in the form 

4 

Ck ^ k = 0, (3) 

k = i 



where, c*, are unknown complex constants. The above equation can be rewritten as 

4 </>* 
c* 

k=i 



E Vk - 



(4) 



The coefficients are determined by choosing four linearly independent plane waves (i.e., 
trial solutions of (2)) traveling along the lines joining the central node to the four neigh- 
boring nodes; i.e., we choose 



'th. — p-jdk cos(a* -a,) 

V>0 



/=!,.. .,4 



( 5 ) 
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where (dk^k) are the polar coordinates of the neighboring nodes with respect to the 
central node 0. The linear system of equations resulting from (4) can then be solved 
for the unknowns c*. It is interesting to note that the values for the coefficients reduce 
to the standard values when the mesh conforms to a standard coordinate system. The 
computational domain is open ended at the sides and one cannot choose a five-point star 
mesh as above. At the two ends of the domain, a modified mesh as shown in Figs. 3a and 
3c is used. For the mesh shown in Fig. 3a, it is possible that the plane waves traveling 
along directions 3 and 4 become linearly dependent, or nearly so. To avoid this degeneracy, 
we reverse the direction of the plane wave traveling along 0 — 3 and solve the linear system. 
Likewise, to avoid degeneracy with the mesh of Fig. 3c, we reverse the direction of the 
plane wave traveling along 0 — 4 and solve for the coefficients. 

At the lower (object) boundary of the computational domain, the impedance boundary 
condition (1) is applicable to the total field H z — x + Vb where x = ^ 0z is the incident 
field. For small distances, the solution of (1) along the normal results in 

4'S = ~xS + « + k " = -xS + « + x")0n, ( 0 ) 

where h n = x - xj) 2 + (yj* - y 0 n ) 2 is the normal distance between scatterer and the 
first interior layer at the nth node, A" is the value of A 5 at the nth node, 
and Xm = XO 

Additional conditions are needed at the artificial, terminating boundary. To simulate 
free-space, boundary conditions proposed by Mei et a/.,[4] are assumed to hold on the 
boundary. Accordingly, on the five-point star mesh of Fig. 3b formed by layer numbers 
M — 1, M, and M + 1, the scattered field satisfies (in the local coordinates) 

4 

tpo + '^2 a ktpk = o, ( 7 ) 

*=i 

where the complex coefficients ak axe different from c*. They are determined in a fashion 
similar to the latter except that the trial solutions are no longer plane waves but are 
dependent on the geometry of the scatterer and the location of the mesh. Such trial 
solutions can be obtained from the near-zone behavior of the field and further discussed 
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in section lid. It is important to note that a discrete equation of the form (7) over a five- 
point computational molecule of Fig. 3 can only capture spatial derivatives upto second 
order (without the cross terms). Equations (3), (6) and (7) are combined to result in a 
matrix equation of the form 

A]^ 3 + + C 1 M' 2 = F 1 (8) 

A n ^ n_1 +B n 'F n + C n ^ n+1 = F n , 2 < n < A — 1 (9) 

A/v ^ -1 + B n V n + C N y N ~ 2 = F n , (10) 



where \F n = [?/>", ^2 >• • • is the column vector of unknown values on the M interior 

layers at the node n, A n ,Bn,C n are banded matrices of order M x M given by 



, 2,n 

0 



A n = 



0 

o 3,n 



0 0 



0 
0 

n M- i,n n M-\,n 



a 2 c 3 



~ «3 1 C 2 



B n = 



1 + P n c\’ n C 2 ’ n 0 



o 3 ,n 



3,n 



c n = 



0 

c 2 4 ' n 0 

0 cl' n 



0 

0 



n Af-i,n n M -\,n n M - 1 ,n 



a 2 c l 



— a\c 2 



a 2 c 2 



0 

0 



n ^M-\,n _ n Af-i,n 
a 3 c 4 



a 2 c 4 



0 0 . 

and F n is the excitation column vector given by 

( -c 2 ,’"(W - xS ) 1 
0 

F n = 

0 



( 11 ) 



( 12 ) 



(13) 



(14) 



In the above matrices, c™’ n is the kth finite difference coefficient associated with the node 
at (a;^,y^), and a% is the Ath MEI coefficient associated with the node at (^”,,2/^)- The 
system of equations defined by (8-10) can be solved efficiently by the Ricatti block-by- 
block elimination technique [8] described in He. After the scattered field is solved in the 
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interior domain, the total field on the scatterer can be recovered from (6). The field at any 
point can then be determined from the total field on the scatterer by using the integral 
representation given in the next section. In the next few sections, we give details of the 
various steps presented in this section. 

lib. Integral Representation: 

Fig. 1. shows an electric line source J 0 located at 0(x 0 ,y 0 ). It is required to find 
the total field (£, i/), at an arbitrary point P (x p ,y p ) over the irregular, inhomogeneous 
terrain. We will employ the reciprocity theorem [2] to express the total field in terms of 
the scattered field and an incident field, (£ 0 ,77 0 ), defined to exist over a flat, constant 
impedance plane, Co, of normalized impedance A 0 . The general relationship between a 
surface impedance A and the local ground constants (e 0 e r , /z 0 , a) is 

A = 1 /\/( e r - jcrr), (15) 

where a T is the relative conductivity representing <r/W 0 , a being the actual conductivity 
of earth. The above relationship is used to evaluate not only A 0 , but also A 5 that will 
be encountered below. The field (E^H) is different from (£o,7/o) due to (a) the terrain 
inhomogeneities, or (b) terrain irregularities, or (c) both. Let be the fields due to 

sources ( J \ , M \ ) located at P(x p , y p ) over the reference impedance plane. The fields ‘O’ and 
‘1’ satisfy the impedance boundary condition y x Eo = rj 0 A 0 y x (y x Ho) on Co, where 
y is the unit normal on C 0 . The total fields satisfy the impedance boundary condition 
y x E = 7] 0 A s h x (n x H ) on the terrain surface C, where the impedance A 5 is possibly 
different from A 0 , and may vary from point to point. The surface C is assumed to consist 
of flat portions C j where the impedance A s deviates from A 0 and/or C t where the terrain 
deviates from being flat. The remaining portion C r of the surface is concident with the 
corresponding portions of Co- The surface Cb is the projection of Ct onto Co- Let Coo 
be a semicircular cylinder bounded at infinity, and A be the interior region bounded by 
the surfaces C and Coo- The unit normal h points into the region A. It is assumed that 
the ground deviates from the reference plane only in the positive y direction. Thus, we 
consider only hilly obstacles and do not address the problem of trenches in a flat plane. 
The ground constants are, however, allowed to have variations along the terrain. 
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The fields (E — Eq, H — Hq ) have no sources in the region A, and satisfy source-free 
Maxwell’s equations, while the fields (Ei,Hi) satisfy Maxwell’s equations with sources 
From a two dimensional version of the reciprocity theorem [2], we have 



/ 



(E-E 0 ) xH 1 -E 1 x (H-H 0 ) 



h dt 



= // 



(E-E 0 )-J 1 -(H-H 0 )-M 1 



da 



C r +C j -\~C t +C oo A 

(16) 

The contour integral over vanishes due to the Sommerfeld radiation conditions [11]. 
Since (E,H) and (E 0 ,Ho) satify the same boundary conditions on C rj the integral over it 
vanishes. Applying reciprocity theorem to the source- free region bounded by C\ — Cfe, we 
obtain 



J [£ 0 x Hi - Ei 

c,-c b 



xH o 



• hdi = 0, 



(17) 



Hence 



J Eq X Hi - El X Hq 
c t 




x Hi - El 



x Ho 



hd£ = 0. 



(18) 



The last equality follows from the fact that the ‘O’ and ‘1’ fields satisfy the same impedance 
boundary conditions on Cb ■ Making use of these in (16) we get 




- E q ) • Ji -(H- H 0 ) ■ Mi 




Cf+C t 



— Ei x H) • h dt. 



(19) 



Next, we pick a i-directed magnetic line source of unit voltage for Mi and choose Ji = 0 
to arrive at 



z ■ H(p) = z ■ H 0 (p 



>-/ 

Cf + C t 



(h x E)- h[ z) + (h x H)- E\ z> dt, 



U) 



( 20 ) 



where (e[ z \ H[^) is the field due to a l-directed magnetic line source of unit voltage 
located at the field point P. We relate (n x H) and (h x E) on Cj + C t by means of the 
impedance boundary condition to finally obtain 

z ■ H(P) = z ■ Ho(P) ~ J [ E\ z \q ) - rioA s (n x ^ (z) (Q))] • (n x H(q)) dt Q . (21) 

C/+C, 
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The above equation is rewritten as 

Hz(p) = H 0 z (p) - f H z (Q)Q(Q,p)d£ Q , (22) 

C=Cj -h C t 

where 

6(q,p) = e[ z ( \q) - r, 0 A ( 23 ) 

is defined as the Green’s function at Q due to a ^-directed magnetic line source located at 
P over a plane of impedance A 0 . The unit tangent t to the contour is defined such that 
£ x h = z. Equation (22) may be used to set up an integral equation for the unknown 
H z and solved numerically. However this will be computationally intensive as already 
discussed and will not be considered in the present report. In the next section, we present 
expressions for H o z and Q(q,p) which allow rapid numerical computation. 

lie. Incident Magnetic Field and Green's Function : 

Fig. 4. shows a transversely polarized electric line source located at (x 0 .y 0 ) over 
a plane of constant impedance A 0 . Due to the impedance boundary condition we have 
Eox = 7j 0 A 0 Ho z . The line source carries a total current I Q and has its current moment 
directed along the unit vector i which makes an angle 9 0 with the x-axis. The current 
moment corresponding to the mirror image of the source about a perfectly conducting 
plane at y = 0 points in a direction i f which makes an angle i x — 9$ with the rr-axis. The 
current density of the line source is expressed as 

J = iI 0 k 2 o S(x - x 0 )6(y - y 0 ) ^ k 2 J i (24) 



where S(-) represents a unit delta function. The magnetic vector potential A = xA x AyA y 
satisfies 

V 2 A + A=-yJ. (25) 



It is easily seen from Maxwell’s equations that 



U3 

h 02 = - 

n 0 



dA y dAj 



E 0x = -jr) 0 



dx dy 
dHoz 
dy 



(26) 

(27) 
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Equation (25) can be seperated into its respective cartesian components. A vertical current 
excites a vertical component, of the vector potential, whereas, a horizontal component 
excites only a horizontal component, A x . It is to be noted here that a horizontally polarized 
line source will only need A x for the satisfaction of the boundary conditions at y = 0. This 
is in contrast to the three-dimensional case where a HED requires both A x and A z for the 
satisfaction of boundary conditions. We take a Fourier transform on both sides of (25) 
with respect to the x-axis, and denote the transformed quantities with a tilde and the 
normalized transform variable by k x . When this is done, an inhomogeneous harmonic 
equation with respect to y is obtained in the transformed domain. The components of the 
vector potential can be determined from it as 



A x (k x , y] x 0 ,y 0 ) = //q/ ° C 2j°/ klZ ° {e~ my - yA + Rnifte-Wy+yA} 
Holo sin 6 0 e jkxX ° 



A y (k x ,y-x 0 ,y 0 ) = 



47rj/? 



je-jOls/— ».l + , 



(28) 

(29) 



where Rh(/3 ) and R v (fl) are the reflection coefficients for the horizontal and vertical polar- 
izations respectively, and k\ + fi 2 = 1. Imposition of the impedance boundary condition 
at y = 0 yields 



R„m = -M0) = 1 - (30) 

The components of vector potential in the space domain are obtained by taking the inverse 
Fourier tranformation of (28) and (29). Letting the free-space Green’s function II as 



n( 



oo 

V) Vo) J 



z -j[kx(z-x o )+0\y-y o \ 

7 t(3 



dk x = H^W{x - x o y + (y - y G ) 2 ], (31) 



we arrive at 



A = |fll(x, y; x 0 ,y 0 ) + i' [II (x,y;x 0 ,-y 0 ) - 2A 0 <5II (x, y; x 0 , y 0 )]| , (32) 

where 

°r e -j[/?(2/+2/o)+M*-*o)] 

6U(x,y,x 0 ,y 0 )= J + A ) ^ 

— OO 

The integral defined in (33) is of the so called Sommerfeld type and difficult to evaluate 
directly. To efficiently compute the incident magnetic field and the Green’s function, we 
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will incorporate the ideas contained in [12] and [13]. For small distances, we use a modified 
method of [13]. For large horizontal distances, where the method of [13] fails, we have used 
a modified version of the method in [12]. 

For small distances, we employ the clever device originally used by Chow, et al , [13] 
to calculate the Green’s function for microstrip structures. The key to the procedure is 
to convert (31) to an integral over the complex /9-plane and make approximations to (33). 
Firstly, we note that the imaginary part of (3 must be chosen appropriately (i.e., G(/3) < 0) 
such that integral in (31) converges for all |y — y 0 1 >0. By a change of variable, (31) may 
be converted to an integral in terms of /? as 

r9 x , 2 f e~ 3 ^ y ~ y °\ 

H ( 0 [\/(x - x 0 ) 2 + (y - y 0 ) 2 ] = - / : cos [k x (x - x 0 )]d/3, (34) 

7T J k z 
r, 

where the contour of integration, Tj, is shown in Fig. 5. Similarly, 

2 f e~^ ( ' y+y °' > 

6U(x,y;x 0 ,y 0 ) = - / cos [k z (x - x o )]d0. (35) 

it J k z [p + A 0 ) 
r, 

The contour of integration for <511 may be deformed to a modified contour, T 2 , since no 
singularities are enclosed between it and Ti (see Fig. 5). The contour T 2 consists of two 
parts: a straight line, drawn from (1,0) to (0, — jTq ), for some chosen positive constant 
To, and a portion coincident with that of Ti. Most of the contribution to the integral 
comes from the former as the integrand decays exponentially on the latter for sufficiently 
large y + y Q . The trick is to approximate the spectral function 1 /(/? + A 0 ) over as a 
finite sum of complex exponentials of the form 



1 

W+ A 0 ) 



n t 

« ^2 Ake~ tkd . 

h=i 



(36) 



The complex constants Ak and can be determined by the Prony’s method [14]. Substi- 
tuting this into (35) we have 



<5n(x,y;x 0 ,y 0 




e -j0(y+yo-jt k ) 

k z 



cos[k z (x — a: 0 )] df3 
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"4fc [■E’l( 2 2 — jtkTo) + E\(z 2 — jtkTo )] • 

7T z — ' 



( 37 ) 



fc=l 



The first term of the last equality follows by comparison with (34). The quantity Rk is the 
distance between the observation point (x, y) and the complex image at (x 0 , —y 0 +jtk) given 
by Rk = \J (x — x 0 ) 2 + (y + y 0 ~ jtk ) 2 J with the square root chosen such that > 0? 

and 22 = [(y + yo) + j( x ~ ^o)]To. The remaining integral in (37) over the negative 
imaginary axis has been performed in closed form in terms of exponential integrals, E i(-) 
[21], by using the approximations (3 + A 0 « (3 and k x « —j[3. The quantity ^11 can 
be calculated much faster using (37) rather than (33) or (35). However, when y + y a is 
small and |x — x 0 | is large, the contribution arising from the tail of becomes important, 
and the exponential approximation to the spectral function in (36) requires many terms. 
Using only a few terms leads to large errors in the computation of <511. The accuracy can 
be greatly improved at the cost of computational time by expressing the spectral function 
1/(13 + A 0 ) over r s as a Laplace integral of the form 



1 

(f3 + A 0 ) 



OO 

J e-^e-^dt. 

t = 0 



(38) 



Such a procedure was suggested in [12]. The Laplace transform of the spectral function in 
(38) converges for — 7r/2 — arg(A 0 ) < arg(i) < 7r/2 — arg(A 0 ). Substituting this into (35) 
and changing the order of integration, we arrive at 



OO 

OT(x,y;x 0 ,y 0 )= J e'*’ 1 H™ (D) dt, 
1 = 0 



(39) 



with the complex distance D = yj(x — x 0 ) 2 + (y + y 0 — iO 2 defined such that 5ft(D) > 0 
for representing an outgoing wave and 5( D ) < 0 for the convergence of the integral. Even 
though it is much faster to compute <5n using (39) than using (35), the integration can be 
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slow for certain combinations of A 0 ,(;r — x 0 ), and (y + y 0 ). This is because the integrand 
decays rather slowly on the positive real axis of the f -plane. To further aid fast evaluation 
of the integral (39), we deform the original contour in the complex t— plane to one over 
which the asymptotic form of the Hankel function decays most rapidly (steepest descent 
path). This is obtained by setting 3?[£)(<)] = = 0)] = r 2 on the deformed contour, 

where r 2 = yf(x — x 0 ) 2 + (y + y 0 ) 2 is the distance of the observation point from the mirror 
image. If t = u + jv, the steepest descent path through t = 0 can be obtained as 

u — + [0, oo), V = -(y + y 0 ) + r 2 + (40) 

y u z + 

on which the distance function is 



D = r 2 — ju] 



■u 2 _ + (y + VoY 
u 2 + ri 



(41) 



Evaluation of the integral over the steepest descent path permits rapid computation 
of 6IL. Although this is true uniformly for all x — xq and y + y 0 . it is less efficient than (37) 
for small values of r 2 . This is due to the presence of the Hankel function with complex 
arguments in the integrand of (39), which, one has to evaluate repeatedly to perform the 
integration. Hankel functions need be computed only for a fixed number of arguments in 
(37) in contrast to the several tens of times needed when (39) is used. It is very important 
indeed to reduce the time required to compute the incident field and the Green’s function, 
particularly when one has to evaluate them several thousands of times as in the present 
method. This point will be appreciated shortly when we discuss the MEI method. 

Combining (26) and (32) the incident magnetic field can be obtained as 
— kol 0 



H oz = 



4 j 



sin(0 o - e 1 )H[ 2 \r l ) + sin(# 0 + 



q g 's 

+ 2A o (sin0 o — + cos 0 0 -^)8H(x,y,x o ,y o ) j, 



(42) 



where rj = yf{x — z 0 ) 2 + (y — y 0 ) 2 , and #i(#2) is the positive angle made with the x-axis 
by the straight line joining the observation point and the source (mirror image) point. 

Using a similar analysis, the Green’s function (defined through (23)) at the point 
Q (x q ,y g ) due to a unit voltage, z-directed, magnetic line source located at P(x p ,y p ) can 
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be obtained as 



G(Q,P) = sin(6 - 0 3 )H {2 \r 3 ) + sin(0 - 6 4 )H[ 2 \r 4 ) - j A s [ff< 2) (r 3 ) + H ( 0 2 \r 4 )\ 

6U(x q ,y q -x p ,y p ) /, 



+2A 0 



d d 

jA s + sin#— cos 9 



dx n 



dy q . 



(43) 



where r 3 = y/{x q - x p ) 2 + (y g - y p ) 2 , r 4 = y/(x q - x p ) 2 + (y, + y p ) 2 , 0 3 (0 4 ) is the positive 
angle made with the z-axis by the straight line joining the observation point Q and the 
source point P (mirror image), and 0 is the angle made with the positive x-axis by the 
unit tangent at Q. 

An integral equation may be derived for the unknown, H z , on the scatterer by substi- 
tuting the above Green’s function into (22) and taking the limit as the point P approaches 
the scattering boundary from outside. When this is done, one gets 



H z (r) = 2H 0z (p)-2 f H z (q)Q(q,p) di Q , p,Qec J+ c t , (22') 

Jc=c ] +c l 

where the integral is understood to be of Cauchy’s principal value type. The first term on 
the right can be regarded as the physical optics approximation to the surface field in the 
illuminated region. 

The far-zone (r — » oo) fields can be determined by using the principal asymptotic 
forms for the various Hankel functions. For the quantity (511, far-zone approximation can 
be obtained by deforming the path of the integral in (33) and evaluating by the saddle 
point method. The result is 



<5n(x,y;x 0 ,y 0 ) ~ J -j- 6 

V 7r r 2 (A 0 + 



(44) 



Having provided the integral representation and expressions for the various fields, we 
present in the next section details on the MEI method of terminating the computational 
domain in the finite scheme. 



lid. MEI Method: 

As already remarked, the MEI method allows one to generate near-field conditions for 
the scattered field rp to simulate free space at the terminating boundary. We will describe 
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the procedure to determine the coefficients a* in (7). Mei et al [4] postulated that there 
exists a linear combination of the scattered field over a small, discrete, spatial domain, p*, 
such that 

K 

— 0. (45) 

k = 1 

The coefficients a * are postulated to be dependent upon the location of the field and 
geometry of the scatterer, but, independent of the excitation of the scattered field. The 
last postulate enables one to determine the coefficients a * using a finite set of linearly 
independent scattered fields caused by different excitations. It is not the purpose of the 
present report to test the validity of these postulates. Rather, we employ this method to 

develop a relatively fast numerical method for predicting both the sky wave as well as the 

ground wave for propagation over inhomogeneous and irregular terrain. The starting point 
is an accurate representation of near fields such as equation (22). The scattered field is 
given from (22) as 

= H,(F) - Ho,(p) = -[ H,(p')g(P,p)de'. (46) 

C 



Making use of this in (45) we see that 




(47) 



The specific excitations (which they termed as metrons ) suggested by Mei et al., were 



H 2 = 1, cos(27rs), sin(27rs), cos(47ts), and sin(47rs), 



(48) 



where s is the normalized arc length that varies between 0 and 1 on the scatterer. We will 
label these as sinusoidal metrons. We were not able to consistently generate meaningful 
results using sinusoidal metrons, and thus considered other choices. A closer look at (47) 
provides insight as to why the sinusoidal metrons cannot be expected to work all the time. 
Since the coefficients a * have been postulated to be independent of the excitation (which 
determines H z (p ') on the scatterer), it is implied that the terms within the bracket be 
identically zero for all points p' on the scatterer. The term C/(p',pk) is the field at p' due 
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to a line source located at pk . Hence we require that the field due to a linear combination 
of the these line sources (J\ + 1 in number) vanish pointwise on the scatterer. But this 
would be difficult to accomplish as there are only & few discrete sources (I\ is typically 4-5 
depending on the type of the computational molecule chosen). This is particularly true 
for a large scatterer. A more reasonable criterion is to require that the field be minimum 
in some sense on the scatterer. We choose the coefficients a * so that integrated square 
residual 



R 



-f 



K 



G(p',po) + ^2akG(p',pk) 



k — 1 



d£' 



(49) 



is a minimum on the scatterer. In other words, we determine the coefficients by requiring 
that = 0. This position was first taken by Jetvic and Lee [16], who considered the 
case of perfectly conducting cylinder to demonstrate the success of this approach. This 
criterion results in 



K 



f G(p',Pk)G*(p',Pn)d£' = - J G(p\po)G*(p',Pn)dt', n = 1, 2, . . . , K. (50) 

fc=1 c c 

The coefficients ajt can be obtained by solving the linear system defined by (50). We shall 
label the particular choice of metrons in (50) as the G* metrons. The coefficient matrix 
in (50) is Hermitian symmetric. The number of integrals required per a five-point star 
mesh ( K = 4) is 14. If every node is utilized in the evaluation of an integral, the total 
number of flops required to perform the integration per a five star mesh is 14iV. The 
number of Green’s function evaluations per a five-point star mesh is 51V. The total time 
required to fill the coefficient matrix in (50) for N nodes on the object boundary is then 
14 N 2 tf + 5 N 2 t g = 0(1V 2 ), where tf and t g are the respective times required per flop and 
a single evaluation of the Green’s function. For example, if it takes 1 millisecond for a 
single evaluation of the Green’s function, the total time spent in the evaluation of the 
Green’s functions for N = 1000 would be 83 minutes. It is therefore very important to 
minimize the time required to calculate the Green’s function. The total time needed to 
fill the matrices (11) through (14) is still of order 0(N 2 ). The fill time may be reduced 
approximately by a factor of half by skipping every other point in the evaluation of the 
integrals. The next section deals with the inversion of the block system defined in (8-10). 
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lie. Inversion by Ricatti Transform 

The overall matrix generated by the block system of equations in (8-10) is highly 
sparse and is almost tridiagonal. It may be inverted efficiently by the Ricatti transform 
method [8], which was originally described in [17]. To this end we write 

*"- 1 = R n *" + S„, n = 2,...,A, (51) 

where the transform matrix R n is of order M x M, and the shift vector S n is of order 
M x 1. Substituting into equation (9), we get 

*" = - [A„R n + B„] -1 C„$ n+1 + [A n R n + B,,]- 1 [F n - A n S„] . (52) 

We then observe by comparing with (51) that for n = 2, . . . , A r — 1 

Rn+i = — [A„R„ + B n ] -1 C„, S n+1 = [A„R„ + B n )~ 1 [F n - A„S„]. (53) 

The end equation (8) may be used to determine the transform matrix R 2 and S 2 . 
Substituting (51) into (8) and carrying out some algebraic manipulations, we get 

R 2 = (A 1 C 2 - 1 A 2 - Bj)" 1 (Cj - A 1 C 2 - 1 B 2 ) , (54) 

and 

s 2 = (A 2 - C 2 A“ 1 B 1 )" 1 (F 2 - C 2 A“ 1 F 1 ) . (55) 

The higher order matrices can be determined from the above two by using the recursive 
relations in (53). Finally, by using the other end equation (10) in conjunction with (51), 
we solve for 'F N as 

= [(A/v + C/vR/v-i)R-/v + B/v] 1 [F a — C/vS^v-i — (A^ + C^R^_i)SAf] . (56) 

The remaining vectors T' V_1 , ^f N ~ 2 , . . . , T fl can then be determined by the tranformation 
equations in (51). Each of the matrices A„, C„, S„, and has M non-zero elements 
while B n has 3 M non- zero elements. The matrix R n is dense and has M 2 elements. The 
total number of non-zero elements, N 3 , that need to be stored with this algorithm is then 

N 3 = (N - 1 )(M 2 + M) + 6MN ~ (M + 7 )MN = O(N). (57) 
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For example, with N = 1000 and M = 5, the algorithm needs around 0.96 MB of RAM for 
implementation in double precision. If N is now increased to 10,000, the required memory 
is increased to 9.6 MB. This is in contrast to the integral equation methods, where, a 
scaling in N by a factor of 10 results in 100 fold increase in memory. The inversion time 
of the algorithm is only of order O(N). The overall CPU time of the present method is 
dominated by the matrix fill time discussed in the previous section. 
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III. NUMERICAL RESULTS 



The entire code was developed in double precision in FORTRAN. The linear equations 
generated by (4) and (5) for the determination of the coefficients c * were solved by Gaussian 
elimination. Due to the generation of a highly ill-conditioned matrix, Gaussian elimination 
is, however, not reliable for the inversion of the matrix associated with the coefficients a 
By treating it as a linear least square problem, we have used Singular Value Decompsition 
(SVD) [19] to construct a minimum norm solution to the matrix equation. SVD is also 
appropriate for the case of sinusoidal metrons where we generate more equations that the 
number of unknowns. The effective rank of the matrix is determined by treating as zero 
those singular values which are less that a predetermined number ‘Rcond’ times the largest 
singular value. The condition number, k 2 , of the matrix in the 2-norm is the ratio of the 
highest singular value to the lowest singular value. Of course, if Rcond is set less than the 
infimum of l//-c 2 over all the nodes on the boundary, SVD uses all the singular vectors and 
produces the same result as Gaussian elimination for a square matrix. Integration in (4G) 
was performed by the Simpson’s rule either using all nodes on the scatterer or using every 
other node. The latter reduces the integration time by a factor of half. The geometry of the 
scatterer was specified in a discrete form by the nodes at which the unknowns are defined. 
Interpolation with quadratic elements [20] was used to generate a continuous object. In 
this way the code was capable of handling a rather arbitrary geometry. The various Hankel 
functions with complex arguments in (37) and (39) were generated by implementing the 
algorithm of du Toit [15]. 

To verify calculations by the complex image approach of section lie, we first present 
results on the computation of the incident magnetic field. Fig. 6 shows the magnitude of 
the normalized magnetic field due to a vertically polarized line source on the surface of a 
flat, lossy plane. The ground constants correspond to e r = 15 and a T — 6. The magnetic 
field is normalized to the free-space value, which is true of all the results shown in the 
report. The constant To was chosen to be at least 5 to make the approximation in (37) 
valid. Fourteen complex images were chosen to approximate the spectral function, although 
10 were found to produce identical results. Expression (37) was used for r 2 < 50, while 
(39) was used otherwise. The upper limit for u on the steepest descent path in the integral 
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(39) is determined such that the imaginary part of the complex distance D takes a value 
of 5. The arc length value in Fig. 6 is equal to zero right under the source. Comparision is 
shown with the calculation based on Sommerfeld integral given in [1]. Excellent agreement 
is seen between the two, thus validating the complex image aproach. Several other cases 
have also been successfully tested including the cases of large conductivity such as e r = 10, 
and a r — 200. Where applicable, we used N t = 10 and T 0 = 5 in the numerical results 
shown in this report. 

We next compare the numerical results for propagation over a semicircular boss on 
a perfectly conducting plane (A 0 = 0), for which the exact solution in terms on cylin- 
drical harmonics exists [1]. Results are shown both for the sinusoidal metrons and the 
G* metrons. A vertically polarized source (0 O = 90) was assumed to be located at 
( — 0.75A, 0.1A) near a semicircular boss of radius 0.5A with center at the origin. Six layers 
(M = 4, recalling that M denotes the number of interior layers not counting the object 
and outer boundaries) with an inter-layer spacing, h, of 0.05A were used to discretize the 
computational domain. The distance between the outermost layer and the object bound- 
ary is (M + 1 )h = 0.25A. The outer boundary was discretized into roughly twenty segments 
per wavelength resulting in N = 49. Rcond for SVD was set at l.D-6. The arc length, S', 
takes a value 0 at the left end of the boss and increases in the clockwise direction. Fig. 
7 shows a good agreement for the magnitude of the surface magnetic fields between the 
numerical and the exact results both for the sinusoidal and G* metrons. The supremum, 
(indicated as ‘Con’ in the plots), of the 2-norm condition number of the 4x4 coefficient 
matrix in the case of G* metrons, and the 5x4 matrix in the case of sinusiodal metrons 
over all nodes on the boundary is also indicated on the plot. Fig. 8 shows a comparison 
of the relative pattern of the source in the presence of the boss. Good agreement is seen 
except near the grazing angles in the shadow region where a discrepancy of about 2 dB 
is seen. The accuracy of the solution can be affected by varying Rcond, M, N and h. 
While it is generally true that the accuracy improves when M or N is increased or when 
is h decreased, difficulty to compute the various quantities precisely over a small, discrete, 
spatial domain tends to obscure this. The last two parameters affect the discretization 
errors, while M and h together influence the distance between the object boundary and 
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the outer boundary. All other paramaters being fixed, if M is increased and h is decreased 
while maintaining a constant seperation between the object and outer boundaries, the so- 
lution improves only slightly as indicated in Figs. 9 and 10. In Fig. 9 we use M = 6 and 
h = 0.0357A, while in Fig. 10 we use M = 8 and h = 0.0278A. The slight improvement 
over the results of Fig. 7 is thought to be due to a decreased discretization error achieved 
by the use of smaller h. Compared to Fig. 7 it is also seen that increases in both cases. 
This is due to the larger size of the matrices involved. The effect of integration on the 
numerical solution was also investigated. In one case, integration was performed by using 
all nodes on the boundary, while in the other case every other node was used. The former 
is labeled as full integration, while the latter as half integration. Fig. 11 shows the results 
with full and half integration. The solution is less sensitive to integration for the G* case 
than it is for the sinusoidal metron case. This trend has been noticed for other shapes 
as well. To see the effect of node density on the solution, the number of segments on the 
outer boundary was increased to 29 per wavelength. This was with a view to increase 
the accuracy of the solution. Fig. 12 shows the solution with the higher node density. 
The G* metron solution appears to have improved slightly, whereas the sinusoidal metron 
case deteriorated slightly when compared with the 20 segment /wavelength case of Fig. 9. 
Note the increase £1 for both cases when compared to Fig. 9. This is because the fields 
calculated over a small spatial domain tend to be more linearly dependent than over a 
larger domain. 

Parametric studies were also made for larger cylinders. Fig. 13 shows the results 
for a radius of 5A. With the indicated choice of parameters, the number of nodes on the 
boundary is N = 331. Once again the agreement between numerical and exact results is 
good. Some spurious oscillations are, however, seen in the numerical solution in the last 
part of the shadow region. The total computational time needed on a 80486-50 Personal 
Computer for the calculation of ground wave data at 331 points and the far-zone data at 
181 angles utilizing every other point for integartion was 2| minutes. Almost all of the 
time was required for filling the various matrices, with the inversion taking only a paltry 
1 second. Fig. 14 shows the relative pattern of the source. It is seen that a slightly better 
agreement is obtained with the sinusoidal metrons than with the G* metrons. As in the 
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previous case the accuracy improved only slightly when the number of interior layers is 
increased without changing the distance between the object and outer boundaries. This is 
clearly seen from Figs. 15 and 16 which show results for M = 6 and M = 8 respectively. 
For the sinusoidal case, however, the result with M = 8 is once again less accurate when 
compared with M = 4. Increasing Rcond had no effect on the solution. However, when 
the number of segments on the outer boundary is increased to 34 per wavelength, a more 
accurate solution is obtained as seen in Fig. 17. At this node density the spacing between 
nodes on the object boundary is the same as h . Fig. 18 shows the numerical solution 
with the higher node density of 34 segments per wavelength for M — 8 and h = 0.0278A. 
For the G* metron case the solution remains practically unchanged from that of Fig. 
16, but the result for sinusoidal metron case is more meaningful with the higher node 
density. The effect of moving the terminating boundary further away from the object 
boundary is seen by examining Fig. 19. For the parameters indicated, the outer boundary 
is at a distance of 0.45A from the cylinder. For the G* metron case, the hump in the 
solution around S = 70 is still present, although an improvement is seen at the left end 
of the boss (the illuminated side). Overall, the G* solution agrees better with the exact 
solution when compared with Fig. 13 where the seperation between the boudaries is 0.25A. 
Increasing the node density to 34 per wavelength did not improve the solution any further. 
Further improvement in the solution can be affected by decreasing h which decreases the 
discretization errors. This is clearly seen from Fig. 20 where the solution with h = 0.035A 
compares better than the one with h = 0.05A. The spurious oscillations can be reduced 
by (i) increasing the seperation between the object and the terminating boundaries, and 
(ii) concurrently discarding those singular vectors that are corrupted by roundoff errors. 
The latter is accomplished by decreasing Rcond during the SVD solution of the coefficients 
dk . In general, it has been found that for a given spatial resolution, the coefficient matrix 
becomes more ill-conditioned (as evidenced by a higher condition number) as the seperation 
between the object and outer boundaries gets larger. Consequently, Rcond has to be set 
a lower value in order to achieve a meaningful solution. Fig. 21 shows the solution for a 
separation of 0.6 A and Rcond=l.D-4. The accuracy is seen to be the best of all the results 
shown so far. If however, only an approximate solution is desired, the parameters of Fig. 
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13 are sufficient. The results seem to suggest that the accuracy depends primarily on 
the distance between the object boundary and the outer boundary provided that various 
quantities are calculated precisely. Twenty segments per wavelength, six layers, and an 
interlayer spacing of A/20 seem to produce reasonable results for most cases. If G* metrons 
are used, integration can be done using every other point. 

Next, we consider an obstacle that has both concave and convex potions. A good 
candidate is the case of a gaussian hill for which results are available in the literature [18]. 
Such a hill has also been considered in [1]. The height of the hill, h g (x ), over the y = 0 
plane is given by 

h,(x) = Ac-U-W'. 

We have choosen A = C — 10A/3, B = 0. The obstacle is taken as the portion of the 
gaussian hill where the height exceeds A/100. For the chosen parameters, this occurs for 
|x| ~ 8A. The total arc length of the obstacle is about 18A. The constitutive parameters 
of the earth are taken as e r = 10, a — lOmS/m at a frequency of 1MHz. The latter 
corresponds to o T = ISO. A vertically polarized source is located at (-50A/3, A/100) to 
the left of the hill. Fig. 22 shows the incident magnetic field, H$ z , on the hill normalized 
to the free-space value. The dashed line corresponds to the field over a flat surface. The 
quantity <511 in (42) has been calculated using (39) since r 2 > 50 on the obstacle. As 
expected, the field matches at the two ends of the obstacle with the field over a flat 
surface. In [18] the solution is obtained via Volterra type integral equation starting from 
parabolic approximation of the Helmholtz equation. By its very nature, the method of [18] 
considers one-way propagation and ignores backscattering. The present method, on the 
other hand, makes no such approximations. Fig. 23 compares the numerical solution of the 
present method with that of [18]. Both the G* and the sinusoidal metron solutions have 
been obtained with 20 segments per wavelength node density, six layers, h = 0.05A, and 
Rcond=l.D-6. For these values of parameters, the number of unknowns on the obstacle is 
361. Full integration was employed in both cases. The solution obtained with G * metrons 
agrees fairly well the results of [18], whereas the one with sinusoidal metrons is highly 
erroneous. The increased field strength on the illuminated side is due to focussing by 
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the concave portion of the hill. The G* metron solution correctly predicts it. No such 
focussing is present in the incident field as can be observed from Fig. 22. Thus, physical 
optics cannot be expected to predict the correct field on the surface. It is to be noted 
that the erroneous result with sinusoidal metrons is not due to round off errors, for the 
G* solution has a higher condition number and still produces a good result. This was also 
checked by repeating the numerical solution with various values of Rcond. The sinusoidal 
metron solution remained erroneous irrespective of the value of Rcond. In an attempt to 
improve the solution, a higher node density of 30 segments per wavelength was also tried. 
However, the solution only worsened. It is concluded that the sinusoidal metrons are not 
satisfactory for arbitrary shapes and that they do not lead to the correct values of the MEI 
coefficients. Choice of metrons is not arbitrary as the originators of the MEI method claim. 
The G* metrons have a physical significance in that they represent the fields generated by 
line sources in an environment compatible with the original problem. No such statement 
can be mode concerning the sinusoidal metrons.- 

It is interesting to compare the merits of the present method with G* metrons vis — 
a — vis the method of [1] which solves the problem using a magnetic field integral equation 
and boundary element method. It took 5 minutes and 20 seconds on a 80486-50 PC to 
produce the results for the ground wave data shown in Fig. 5 as well as to generate sky 
wave data at 181 angles. Identical results were obtained with half integration which took 
only 3 minutes. This is in contrast to the method of [1] which took 75 minutes for the same 
problem. Of course [1] uses the free-space Green’s function in contrast to the half-space 
Green’s function employed here. The number of unknowns in [1] was 725 (obtained with 
12 nodes per wavelength). For the same number of unknowns, the present method, being 
of order 0(N 2 ), would take about 21 minutes, which is still faster by a factor of about 
four. If half integration is used the method is faster by a factor of eight. The savings in 
memory in the present method are tremendous. For the problem at hand, the method of 
[1] with 725 unknowns requires a storage space of at least 8.4 MB compared to only 254 
kB needed with the present method. To speed up the calculation of the Green’s function 
we have assumed A„ = 0 in (43). This amounts to assuming that the terrain outside the 
obstacle is made up of a perfect conductor. This results in substantial savings in time. 
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For example, it took 3 hours and 10 minutes on a 80486-50 PC to solve the problem of 
gaussian hill using the exact Green’s function and half integration, compared to only 3 
minutes with the approximate Green’s function! The results were identical except near 
the right end of the hill where the two differed slightly. The value at x — x Q = 160 was 
1.2 with the exact Green’s function compared to about 0.9 with the approximate Green's 
function. One would expect the use of approximate Green’s function to be the worst when 
A 0 is significantly different from zero. We considered the case of highly lossy earth with 
e r = 5 and a = 0.056mS/m at 1 MHz (a r — 1). A gaussian hill with A = \,B = 0, and 
C = 0.76A was considered. A vertically polarized source is placed at (— 3A,A/10). Six 
layers were chosen with h = 0.05A and integartion was performed utilyzing all the nodes. 
Fig. 24 shows the calculations with the exact and the approximate Green’s functions. The 
field that would exist on flat earth is also shown. Calculations with the exact Green's 
function took 30 seconds while those with the exact Green’s function took 45 minutes. 
However, it is seen that the results do not differ much from each other, justifying the use 
of the approximate Green’s function. 

It is possible to increase the accuracy of the numerical results shown in Fig. 23 by 
increasing the seperation between the obstacle and the outer boundary. Fig. 25 shows 
calculations with 6, 8, and 10 layers all with h = 0.05A. It is seen that the result with 
10 layers is generally in best agreement with the results of [18]. The numbers indicated 
within the paranthesis in the caption are the values of Rcond assigned in each case to the 
first 15 wavelengths and the last 3 wavelengths respectively on the boundary. As discussed 
previously, the condition number increases generally as the seperation is increased. We 
selectively filter out the singular vectors which are highly corrupted by roundoff errors by 
choosing a higher value of Rcond. The far-zone fields are however insensitive to the slight 
variations of the surface fields. Fig. 26 shows the relative pattern of the source with 6 and 
8 layers. It is seen that the two cases produce virtually the same results. The radiation 
pattern of the source in the direction of the hill is highly perturbed by the presence of 
the latter. Away from the hill the relative pattern in the presence of the hill is not very 
different from the pattern of a vertical source in free-space except near grazing angles 
where it has a deep minimum. 
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Finally, we present results on the effect of relative location of the source with respect 
to the obstacle on the fields. A vertically polarized source was placed near a gaussian hill 
having A = 1.2A ,B = 0, and C — A. The total arc length of the obstacle is 5.1A. The 
ground constants were e r — 14.9, a r — 14.58. Six layers with 20 segments per wavelength, 
and h = 0.05A were chosen in the numerical model. In one case the source was placed at 
the bottom of the hill at a horizontal distance of 4.53A from the peak. In the second case 
the source was placed at the top of the hill ( x 0 = 0). Fig. 27 shows the normalized ground 
wave on the hill. Once again there is focussing on the illuminated side of the hill when 
the source is at the bottom. The field for the source at the top of the hill is symmetric 
as expected. Note that the field in the shadow region of the hill for xo = — 4.53A is as 
strong as the corresponding field for x 0 = 0. There is no apparent advantage of siting the 
antenna at the top of the hill to receive the ground wave. However, the radiation pattern 
of the source from <j> = 0 to 90° is significantly affected by the presence of the hill. This 
can be seen from Fig. 28 which shows the relative pattern for the two source locations. 
Irrespective of the source location, the radiation pattern has a deep minimum near grazing 
angles as is characteristic of vertical antennas over lossy earth. 



26 



IV. SUMMARY 



A fast, finite difference method that includes all aspects of wave phenomenon such 
as reflection, refraction, diffraction, and backscattering is presented for predicting two di- 
mensional propagation over inhomogeneous, irregular terrain. The terrain is characterized 
by its elevation and impedance, which, in turn depends on the ground constants of the 
earth. Both of these may vary with distance. The terrain topology data is specified at 
discrete points. Interpolation using quadratic elements is done to define a continuous ge- 
ometry. The computational domain for the problem consists of the area bounded by the 
terrain at the bottom and a truncating boundary at the top. To simulate free-space on the 
truncating boundary, discrete, near-field radiation condition of Mei type, derivable from 
an integral representation of the fields, is imposed. Green’s function for half-space is used 
to reduce the number of unknowns. Unknowns are distributed on the terrain only where 
its elevation is non-zero and/or where its impedance differs from a reference value. The 
computational domain is discretized using interior layers between the truncating and the 
object boundaries. Finite difference coefficients valid for an irregular, non-othogonal mesh 
are presented. Accurate expressions are provided for the Green’s function and incident 
fields over a constant-impedance, flat plane. The expressions permit rapid computation as 
they do not involve the troublesome Sommerfield integrals. Results are presented for the 
ground wave as well as the sky wave. 

The truncating boundary can be in the near-field of the obstacle; as near as a A/4 
away from it, and the method works both for concave and convex geometries. Good results 
have been obtained with a node density of about 20 per wavelength and an interlayer 
spacing of about A/20. The metrons used for the determination of the Mei coefficients are 
proportional to the complex conjugate of the Green’s function and fully accomodate the 
environment of the problem. Singular value decomposition which permits filtering of the 
space spanned by the corrupted singular vectors is used to solve for the coefficients. 

The method is attractive for large terrain obstacles where other methods tend to 
be slow. The overall computational time of the method is of order 0(N 2 ). Storage re- 
quirements are of order 0(N), where N is the total number of unknowns (nodes). As an 
example, for a terrain obstacle extending over 18 wavelengths, the number of nodes on the 
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boundary is 361. To compute ground wave data at all of these nodes as well as to compute 
the far-fields at 181 angles, the method takes around 3 minutes on a 80486-50 PC. The 
total storage required for the problem is around 254 kB. It would take around 300 minutes 
with a storage space of around 2.5 MB to solve the problem of a terrain where the obstacle 
extends over 180 wavelengths. At a frequency of 1 MHz the terrain obstacle corresponds 
to one whose arc length is 54 km. At a frequency of 1 GHz it corresponds to an arc length 
of 54 m. 

Certain spurious oscillations and over estimation of the field have been observed in 
some regions of the shadow region. This is attributed to the fact that it is difficult to 
calculate near-fields precisely over an electrically small, discrete domain. These extraneous 
effects can be some what reduced by increasing the seperation between the obstacle and 
the truncating boundary. It is speculated that they can be further reduced by choosing a 
computational molecule that accomodates higher order derivatives. More efficient schemes 
of evaluating the Green’s function and the near-fields will have to be made to further 
reduce the computational time. 

Extension to the three dimensional case is trivial in principle, although not so com- 
putationally. Savings over integral equation methods will be even more dramatic in the 
case of three dimensional obstacles. Although the Green’s function using complex images 
will be little more involved than for the two-dimensional case, it will have the computa- 
tional advantage of being expressible in terms of an exponential functions instead of Hankel 
functions with complex arguments. The fields decay more rapidly with distance in three 
dimensions than in two dimensions. As a result integration could be done little more effi- 
ciently in the former case. The corresponding computational molecule will now consist of 
seven nodes instead of five in two dimensions. Finite difference and terminating boundary 
conditions will have to be developed for vector fields instead of scalar fields. These task 
are presently being pursued. 
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Surface Magnetic Field on Lossy Flat Earth 
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FIG. 6. 
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Surface Magnetic Field on Circular Boss 
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FIG. 13. 
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Source Pattern Near Circular Boss 
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Surface Magnetic Field on Circular Boss 
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Surface Magnetic Field on Circular Boss 
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Surface Magnetic Field on Circular Boss 
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Surface Magnetic Field on Circular Boss 
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Surface Magnetic Field on Circular Boss 
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Surface Magnetic Field on Gaussian Hill 
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FIG. 23. 
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Surface Magnetic Field on Gaussian Hill 
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Pattern of Vertical Source Near Gaussian Hill 
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